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PRIORITY CLAIM 

The present application claims priority to U.S. Provisional Patent Application 
10 Serial No. 60/392,953, filed July 1 , 2002 and U.S. Provisional Patent Application Serial 
No. 60/426,602, filed November 1 5, 2002. The disclosure of both of these documents is 
incorporated herein by reference. 

TECHNICAL FIELD 

15 The subject invention relates generally to optical methods for inspecting and 

analyzing semiconductor wafers and other samples. In particular, the subject invention 
relates to improvements in analysis of the measured optical signal characteristics from a 
sample to determine parameter values for that sample. 

20 BACKGROUND OF THE INVENTION 

As geometries continue to shrink, manufacturers have increasingly turned to 
optical techniques to perform non-destructive inspection and analysis of semiconductor 
wafers. Techniques of this type, known generally as optical metrology, operate by 
illuminating a sample with electromagnetic radiation and then detecting and analyzing 

25 the reflected energy. Scatterometry is a specific type of optical metrology that is used 
when the structural geometry of a subject creates diffraction (optical scattering) of the 
incoming probe beam. Scatterometry systems analyze diffraction to deduce details of the 
structures that cause the diffraction to occur. In general, optical metrology measures 
physical dimensions and/or material properties of structures on a wafer from the optical 

30 response of the wafer, that is, how the wafer changes light that illuminates it. 

Various optical techniques have been used to perform optical scatterometry. 
These include broadband spectroscopy (U.S. Patent Nos. 5,607,800; 5,867,276 and 
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5,963,329), spectral ellipsometry (U.S. Patent No. 5,739,909) single-wavelength optical 
scattering (ILS. Patent No. 5,889,593), and spectral and single-wavelength beam profile 
reflectance and beam profile ellipsometry (U.S. Patent No. 6,429,943). Scatterometry, in 
these cases generally refers to optical responses in the form of diffraction orders 
5 produced by period structures, that is, gratings on the wafer. In addition it may be 

possible to employ any of these measurement technologies, e.g., single-wavelength laser 
BPR or BPE, to obtain critical dimension (CD) measurements on non periodic structures, 
such as isolated lines or isolated vias and mesas. The above cited patents and patent 
applications, along with PCT Application WO03/009063, US Application 2002/0158193, 
10 US Application 10/243245, filed Sept 13, 2002, US Application 2001/0051856 Al, PCT 
Application WO 01/55669 and PCT Application WO 01/97280 are all incorporated 
herein by reference. 

Metrology uses three classes of parameters and their relationships. Structures on 
the wafer have physical parameters, e.g., the thickness of a layer, the widths of a line 

15 structure at various heights (measured generally perpendicular to a face of the wafer, or 
the complex optical index of a material. Most scatterometry measurements are 
performed over a range of independent parameters. Examples independent parameters 
are wavelength (for spectroscopic systems) or angle of propagation. Much of the 
following is written as if wavelength is the independent parameter. The goal of 

20 metrology is to relate measurements to a model of what is on the wafer, where the model 
has parameters that in some sense reflect the physical parameters on the wafer. In some 
cases, model parameters represent exactly the physical parameters in question, or they 
may be related through some mathematical transformation, e.g., the physical widths of 
adjacent periodic lines and spaces may be modeled as a period and ratio. 

25 Most scatterometry systems use a modeling approach to transform empirical data 

into tangible results. For this type of approach, a theoretical model is defined for each 
subject that will be analyzed. The theoretical model predicts the output (reflected) 
electromagnetic field that is generated when an incident field is applied to the subject. 
The theoretical model is parameterized and each parameter corresponds to a physical 

30 characteristic of the subject such as line width or layer thickness. A profile is a collection 
of line widths at various heights, and is thus a collection of parameters. A regression is 
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performed in which the parameters are repeatedly perturbed and the model is repeatedly 
evaluated to minimize the differences between the modeled results and results that are 
empirically obtained. The differences are typically calculated over the range of the 
independent parameters, and then an average difference, such as a squared sum, is 
5 calculated as a single difference. Various norms or other techniques are suitable for 
collapsing the multiple differences into a single working difference. When, the 
minimization reaches some stopping criterion, it is assumed that the model and its 
associated parameters accurately reflect the subject being analyzed. One such stopping 
criterion is that the difference reaches some predetermined level, e.g., that a goodness-of- 

10 fit criterion is exceeded. Another criterion is reached when the reduction of the 
difference from repeat to repeat becomes sufficiently small. Others are possible, 

Evaluation of theoretical scatterometry models is a complex task, even for 
relatively simple subject. As subjects become more complex, e.g., having more 
parameters, the calculations can become extremely time-consuming. Even with high- 

15 speed processors, real-time evaluation of these calculations can be difficult . This is 
problematic in semiconductor manufacturing where it is often imperative to detect, 
processes that are not operating correctly. As the semiconductor industry moves towards 
integrated metrology solutions (i.e., where metrology hardware is integrated directly with 
process hardware) the need for rapid evaluation becomes even more acute. 

20 A number of approaches have been developed to overcome the calculation 

bottleneck associated with the analysis of scatterometry results. Many of these 
approaches have involved techniques for improving calculation throughput, such as 
parallel processing techniques. For example, co-pending application Serial No. 
09/906,290 filed July 16, 2001 describes a system in which a master processor distributes 

25 scatterometry calculations among a group of slave processors. This can be done by as a 
function of wavelength, so that each slave processor evaluates the theoretical model for 
selected wavelengths. The other slave processors will carry out the same calculations at 
different wavelengths. Assuming there are five processors (one master and four slaves) 
and fifty wavelengths, each processor will perform ten such calculations per iteration. 

30 Once complete, the master processor combines the separate calculation and performs the 
best fit comparison to the empirical results. Based on this fit, the master processor will 
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modify the parameters of the model (e.g. changing the widths or layer thickness) and 
distribute the calculations for the modified model to the slave processors. This sequence 
is repeated until a good fit is achieved. 

This distributed processing approach can also be used with multiple angle of 
5 incidence information. In this situation, the calculations at each of the different angles of 
incidence can be distributed to the slave processor. Techniques of this type are an 
effective method for reducing the time required for scatterometry calculations. At the 
same time, the speedup provided by parallel processing is strictly dependent on the 
availability (and associated cost) of multiple processors. Amdahl's law also limits the 

10 amount of speedup available by parallel processing since serial portions of the program 
are not improved. At the present time, neither cost nor ultimate speed improvement is a 
serious limitation for parallel processing techniques. As the complexity of the geometry 
increases, however it becomes increasingly possible that computational complexity will 
outstrip the use of parallel techniques alone. 

1 5 Another approach to rapidly evaluating scatterometry measurements is to use pre- 

computed libraries of predicted measurements. This type of approach is discussed in 
PCT application WO 99/45340, published September 10, 1999 as well as the references 
cited therein. In this approach, a library of expected results is constructed by repeatedly 
evaluating the theoretical model for range of different parameters. When empirical 

20 measurements are obtained, the library is searched to find the best fit. 

The use of libraries speeds the analysis process by allowing theoretical results to 
be computed once and reused many times. Of course, libraries are necessarily limited in 
their resolution and can contain only a finite number of theoretical results. This means 
that there are many cases where empirical measurements do not have exact library 

25 matches. In these cases, the use of a library represents an undesirable choice between 
speed and computational accuracy. 

To overcome this limitation, U.S. Patent Application Publication No. 
2002003 8 196 A 1 (incorporated in this document by reference) describes a database 
method of analysis of empirical scatterometry measurements. The database method is 

30 similar to the library approach in that it relies on a stored set of pre-computed 
"reflectance characteristics." (The reflectance characteristics could be simulated 
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reflectance signal spectra, although in the preferred embodiment they are complex 
reflectance coefficients from which the reflectance signal can be computed without 
resorting to complex electromagnetic simulations). In this case, however an interpolation 
method is used in combination with the database- stored characteristics, making it 
5 possible to achieve measurement resolution and accuracy much better than the database 
sampling resolution. Both the database size and computation time are consequently 
greatly reduced relative to library-based methods. 

A critical element of the database interpolation method is the interpolation 
algorithm itself. In the application just cited (i.e., U.S. Patent Application Publication 

10 No. 20020038196A1), two preferred algorithms are disclosed: multilinear and 

multicubic. Multilinear interpolation is very fast, but has poor accuracy. Multicubic 
interpolation is much more accurate, but can be slow, especially when many parameters 
are being simultaneously measured. In practice, selection of between multilinear and 
multicubic method is done on a pragmatic basis depending on the degree of accuracy and 

15 speed required. Often, this choice is totally acceptable. This is true, for example, when 
the number of parameters is relatively small. Still, it is clear that interpolation methods 
that provide increased speed and accuracy would be an improvement over current 
techniques. 

20 SUMMARY 

An embodiment of the present invention provides an interpolation method for 
analyzing empirical scatterometry measurements. The interpolation method is a variant 
type of cubic interpolation algorithm that exhibits very good accuracy-— comparable to 
multicubic-— but which has much better runtime performance (e.g. for one particular 
25 example the runtime for a six-parameter measurement was reduced from 6.5 sec to 1 .9 
sec). 

The objective of the interpolation algorithm is to estimate values of a reflectance 
characteristic function: 

30 wherein the x's are measurement parameters (e.g., film thicknesses, grating 

linewidths, etc.). The estimate is obtained from a set of pre-computed function values 
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sampled on a rectangular grid in an N-dimensional "parameter space." The reflectance 
characteristic function /[x 0 , x[ , . . . , Jc Ar _ 1 ] is approximated as a linear combination of 

function values selected from the precomputed set, wherein the specific set of selected 
values and associated linear combination coefficients are determined from the parameter 

values x 0 ,x 1? ... The interpolation method requires (N + 1)2 N sample points per 
interpolation (e.g.: 448 for six parameters). This is comparable to multilinear 
interpolation which uses 2 N sample points (e.g.: 64 sample points for six parameters) 
and much better than multicubic interpolation which requires 4 N sample points (e.g.: 
4096 sample points for six parameters). 

The essence of the improved method can be summarized as follows: The ' 
reflectance characteristic function f[x 0 , x, , . . . , x^, ] is approximated as a polynomial 

function of x 0 , x 1 , . . . consisting of a linear combination of elemental monomial functions. 

The monomial terms include all of the multilinear monomials, 

Po P) 

0 1 

wherein each of the exponents p 0 . is either 0 or 1 . (There are 2 N such 

terms.) The monomials also include the multilinear terms multiplied by the square of any 
one of the parameters, 

(x 0 0 x, 1 . . .)xq 

Including all such terms, the total number of monomials making up the 
polynomial interpolation function is: (N + 1)2*. The same number of conditions is 
required to determine the polynomial coefficients. These conditions are determined by 
specifying the function values and its gradient at the corners of an N-dimensional box (a 
database "grid cell") enclosing the interpolation point {x 0 *x,,,..} . (If the interpolation 

point is outside of the database sampling range there will be no grid cell enclosing the 
interpolation point. In this case the grid cell nearest the interpolation point may be used 
to extrapolate the function outside the database range.) The gradient values may be 
obtained by either storing gradient information in the database or by applying a finite- 
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difference derivative estimator to the database-sampled function values. In either case, 
the number of database quantities (function values and possible derivative values) 
required to determine the polynomial coefficients is (N + 1) 2 N . 
The reduced multicubic interpolation function is: 

JoJi,.-MQA}\ *e{0,l,...,jV-l} J 

0<X.<\ 



m 



With the boundary constraints: 

d 

Jk 

The right-hand sides of the constraint equations are predetermined values. 
10 (Number of unknown a and b coefficients = number of constraints = (N + 1)2* .) 

Implicit solution for the a and b coefficients in terms of interpolation basis functions (U 
mdV): 

/[x 0 ,x p ...] = £ f hh U Jq Jx [x 0? x 1? .,.]+ £ #A;y 0 ,y 1>; .r^, 7j ,.„K^i-.-] 

The basis functions have the same functional form as / and satisfy the defining 
15 conditions. 

For x 0 ,jc,,... = 0orl: 

V h J} t Jx 0 , x, , . . .] . = 0 , except U Jo Ji [j 0 , j\ ,...] = 1 
d 

V kj B j [*o,v- ] = o 

20 ^^oj„..[ x o^P--] = o, except— ^ JoJi> ...[y 0 ,y„...] = i 

The following defining relations for U and V satisfy the above conditions: 

U 00t .[x 0 ,x i ,...] = (\-x 0 )(l-x l )---(l + x 0 (l-2x 0 ) + x l (\ -2x, ) + ...> 
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10 



*W . . [*0 ■ ■ "I = *0 0 - *0 f 0 - X ! ) 0 ~ X 2 ) • 



^A;0,0,...[ JC 0'- X: ]'- ••] — ^0; 0, 0,. . . f^O ' ^1 » " • *] ~ 



Jm ^ 1 Jm ■ 1 A ;;/ . ^ 1 A », 



;7o;7„.i x o^,v.-]| j k ^\-j k .-■ F *;7.,7„...K^p--0| x ^ 1 _ JCa 



Combining the above conditions, the following result is obtained: 

t/ A , 7i jx 0 ,x 1 ,...]=(i-x;)(i-x 1 ')-(i+x;(i-2x;)+x l '(i-2x 1 ')+...) 

and: 

h,x„...]=(-rx;(i-^)(i-,;)(i-x;)(i-,;)... 

15 wherein: 



x' = 



As noted above, the interpolation method dramatically reduces the number of 
values that must be used during the interpolation step associated with computing the 
reflectance characteristics of a sample. For example, the number of values necessary to 

20 perform this interpolation when applied to a database wherein each point is defined by six 
parameters is only about one-tenth of the number of values needed with the multicubic 
approach. This computational savings increases for higher parameter problems (e.g.: 
with a seven parameter problem, the number of values needed for the interpolation is 
only one-sixteenth of the values needed for the multicubic approach). Although number 

25 of values used for the interpolation 's reduced, the interpolation method nonetheless 
produces interpolation results that fit all data points, with a smooth curve between data 
points and a smooth transition around data points. In addition to analyzing measured 
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data, the interpolation method can also be used to reduce the complexity and time needed 
to generate a database. 

BRIEF DESCRIPTION OF THE DRAWINGS 

5 Figure 1 is a diagram of a prior art optical metrology system. 

Figure 2 is a perspective drawing of a typical scatterometry subject. 
Figure 3 is a flowchart showing the steps associated with a representative 
embodiment of the interpolation method provided by the present invention. 

1 0 DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 

The present invention provides several related reduced multicubic interpolation 
methods. Each of these methods is specifically intended to be used during the analysis of 
empirical scatterometry measurements. At the same time, it should be appreciated that a 
diverse range of other applications are also possible, including digital terrain mapping 

1 5 and color interpolation for computer graphics. 

Figure 1 shows a scatterometer of a type suitable for use with the reduced 
multicubic interpolation method. As shown in Figure 1, scatterometer 100 includes an 
illumination source 102 that creates a mono or polychromatic probe beam. The probe 
beam is preferably focused by one or more lenses 104 to create an illumination spot on 

20 the surface of the subject under test 1 06. A second lens (or lenses) 1 08 transports the 

diffracted probe beam to a detector 110. Detector 110 captures the diffracted energy and 
produces corresponding output signals. A processor 112 analyzes the signals generated 
by detector 1 10 using a database 1 14. 

As shown in Figure 2, scatterometry subject 106 typically includes a scattering 

25 structure 202 formed on a substrate 204. For the specific example of Figure 2, the 

scattering structure is a grating composed of a series of individual lines. In general, the 
scattering structure may be periodic (as in the case of Figure 2) or isolated. Isolated 
structures include, for example individual lines, or individual vias. The scattering 
structure of Figure 2 is uniform (i.e., exhibits translational symmetry) along the Y axis. 

30 For this reason, this particular scattering structure is considered to be two-dimensional. 
Three dimensional scattering structures are also possible both in isolation (e.g., a single 
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via) or periodically (e.g., a pattern of vias). The seatting structure is covered by an 
incident medium that is typically air but may be vacuum, gas, liquid, or solid (such as an 
overlaying layer or layers). One or more layers may be positioned between the scattering 
structure and the substrate. Alternatively, more than one layer may have a scattering 
5 structure, as is described in U.S. Patent Application Publication Number: 20020158193 
(incorporated in this document by reference). During analysis, a probe beam is directed 
at the scattering structure. For most applications, the probe beam intersects the scattering 
structure at a normal azimuthal angle — the propagation vector is in the plane 
perpendicular to the lines from which the scattering structure is formed. It is also 

10 possible to use a non-normal azimuthal angle. This is referred to as conical scattering. 

The relationship between the incoming and diffracted probe beam is controlled by 
the optical response characteristics of subject 106. To analyze the empirical 
measurements made by scatterometer 100 (i.e., the empirical optical response 
characteristics) a modeling process is typically used. The modeling process is based on a 

1 5 parametric model of the particular subject being analyzed (e.g, Lifeng Li, "A modal 

analysis of lamellar diffraction gratings in conical mountings", Journal of Modern Optics, 
1993, Vol. 40, No. 4, 553-573). The model is evaluated to compute predicted or 
theoretical optical response characteristics for the subject. The theoretical and empirical 
optical response characteristics are compared to determine if the model matches the 

20 empirical results. The model is then perturbed (i.e., the parameters are varied) and re- 
evaluated until the theoretical results and empirical results match within a desired 
goodness of fit. At that point, the parametric model is assumed to match the subject 
being analyzed. 

Database 1 14 is constructed to include a series of pre-computed theoretical optical 
25 response characteristics. Each of these theoretical optical response characteristics is 
stored in association with its corresponding modeling parameters to define an 
interpolation point within database 114. Typically, each stored optical response 
characteristic includes a series of complex reflectance coefficients or scattering matrices 
associated with different illumination wavelengths, incidence directions, and polarization 
30 states, but all associated with the same diffractive structure configuration (materials and 
geometry). The computational representation of the associated structure geometry (e.g., 
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profile shape) is not required for subsequent measurement processes and need not be 
stored in the database. 

Database 1 14 is necessarily limited to include a finite number of interpolation 
points. The reduced multicubic interpolation methods provided by the present invention 
5 provide a mechanism for accurately interpolating optical response characteristics that fall 
between interpolation points within database 1 14 (or in some cases are outside of the 
limits of database 114). This makes database 114 appear to have an infinite number of 
interpolation points without requiring that the theoretical model be reevaluated for each 
optical response characteristic that is not specifically associated with an existing 

10 interpolation point. 

In this document, interpolation (including the reduced multicubic interpolation 
methods provided by the present invention) is described in terms of a function Y[X] . For 
scatterometry, Y[X] is the function that is used to generate theoretical optical response 
characteristics. The interpolation points in database 1 14 are evaluations of Y[X] stored 

15 in combination with the corresponding values of X. The output of Y[X] is typically a 
vector or matrix that typically includes (e.g., for the case of scatterometry) computed 
reflectance quantities such as complex reflectance corresponding to multiple 
wavelengths, incidence angles, and polarization states. The input to Y[X] is typically an 

N -dimensional, vector argument X = (X V X 2 ,...X N ) (For clarity square braces [...] are 
20 used as function argument delimiters and round braces (...) are reserved for grouping; 
vectors and matrices are represented in bold type). 

In the preferred embodiment, Y[X] is sampled at the points of a rectangular grid 

on which each X } takes on the discrete values 

X j =Xf\Xf,...xf' ) ; j = \...N (1) 

25 (The X j sampling values are sorted in increasing order: Xf ] < Xj ] <...xf* j) .) 

Some interpolation algorithms require a specification of not just Y , but also some 
set of partial derivatives of Y at the grid points. In this document, the j-th partial 
derivative operator is denoted 3 . , and product notation is used to denote derivative 

operator composition as follows: 
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d k Y[(X t ,X 2 ,...X N )} = ^Y[(X ] ,X 2 ,...X N )] 

j 

( \ 

Y\dj k Y^d l k, d 2 k2 ...Y 



(2) 



(Note: for A; = 0, d k Y = Y). 
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The grid sampling-point vectors are denoted as: 

X {m) ={X[ nh \X[ m i\...Xp ) )\ mj=Q..M/ 9 j = l...N (3) 
wherein m is a "grid index vector", m = (m 1? m 2 . . m N ) . The interpolation domain is 

partitioned into N -dimensional, rectangular "grid cells" with cell vertices at the grid 
5 points. For a grid index vector m , "cell m " denotes a grid cell defined by the 

constraints Xp ] < X } < X { p +X) ( X {m) is the < c base vertex" of cell m ). Within cell m \ 
Y[X] is approximated by an interpolation flinction Yj {m) [X] : 
Y[X]~Y} m) [X]- 

(4) 

X^ } <X.<X^ +1) ; m,=0...A/,-l; j = !.-.. N 

The collection of cell-bounded interpolation functions Yj (m) [X] defines a global 

1 0 interpolation function Y f [X] over the entire interpolation domain ( Xf* < X j < X { ™' } ). 

The computation of Yj m) is simplified by applying a linear transformation to X 
in order to transform grid cell m into a unit cell with its base vertex at the origin. The 
transformed X argument, denoted x == (x, , x 2 , . . . x N ) \ is defined by 

X — Y^ m ^ 

1 5 The functions Y and F/ yw) translate to transformed functions y and j> y , 

j;W = F[X], (6) 
The Y derivatives transform as follows, 

djylx^iXf^-X^djYiX] (7) 

(and similarly for Y } ). 
20 The vertices of cell m are at X = X (m+5) ; wherein 5 is the integer vector 

s = ,^2, . . ,s N ) with element values 0 or 1 : Sj=0 9 \; j = 1 . . . N . The Equation 5 

transformation maps X = AT (mts) onto jc = s ; thus the problem is to find an interpolating 
function y f [x] approximating y[x] within the interpolation cell 0$jc. <1 (j = 1...N) ? 
given a sampling of y and (possibly) some set of its derivatives at the cell vertices 
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x - s ; Sj =0,1. (In some interpolation algorithms, finite-difference derivatives are 
employed using some external y points sampled outside of the interpolation cell. In this 
case it may be necessary to "pad" the data sampling range with points outside the 
interpolation range to cover these external points.) 
5 Multilinear Interpolation 

This section describes a prior art multilinear interpolation technique. This 
technique is presented because it provides a relatively simple example of the 
nomenclature and formalisms just described. The description of multilinear interpolation 
is also important in later sections which evaluate the performance of the interpolation 
10 methods provided by the present invention. Multilinear interpolation provides an 

interpolation function y } [s] that matches the function being interpolated at each grid 
point of s : 

yi [s} = y[sl 5, =0,1; ..N (8) 

The interpolation function has the form: 
15 y,[x]= X u s[x]y[s] (9) 

■v.s. -0.1 ( /-I...A') 

The Uj[x] factors are "interpolation basis functions" defined by: 

no-*,) (io) 

7 = 1... N. 

in which. x'j is shorthand for: 



f x i9 s. =0 



20 
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Multicubic Interpolation 

This section describes the prior art multicubic interpolation technique. This 
technique is presented (once again) as a relatively simple example of the nomenclature 
and formalisms used throughout this description. The description of multicubic 
interpolation is also important in later sections which evaluate the performance of the 
interpolation methods provided by the present invention. Like the multilinear 
interpolation function, the prior art multicubic interpolation function y } matches y at the 
grid points s . In addition, the first-order derivatives (3, y E , d 2 y } , etc.), and mixed-first- 
order derivatives (d } d 2 y } , d } d 2 d 3 y } , etc.) of y } also match y at the grid points s : 



Y 



( 



n <v riM- ii a/< 



(12) 



kj = 0. 1; *,=0,1; j = \...N 



The interpolation function is 



A;*, =0,1 (j-!... AT) s:s-j=0,i(j=rl...N) 



n v 



(13) 



wherein k = {k x ,k 2 ,...k N ) . The interpolation basis functions, u k s [x], are defined by 



(14) 



15 in which 



4W= 



f(l-x,) 2 (l + 2x,), *,=0 
1(1-*/*,, 



(15) 



and x'j is defined as: 



Xj, 



s r =0 



l-Xj, Sj=\ 
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Reduced-Cubic Interpolation with Exact Derivatives 

As provided by the present invention, reduced-cubic interpolation matches y { and 
its simple derivatives to y at the grid points s (in this document, the term simple 
derivative refers to derivatives that are first order and are not mixed-partial derivatives— 
5 e.g., 9, y } and d 2 y } but not d x d 2 j> 7 or d, 2 y f ): 

yAs]^y[sl yi [s] = d i y[s\; Sj =0 9 l; iJ = \...N (16) 
The interpolation function is: 



10 



15 



( ^ 



i=\...N 



(17) 



The interpolation basis functions, u s [x] and v ( Jx] , are defined by: 



u s [x] = 



^=l..> 7=1. ..AT 

( \ 



(18) 



v,J-v|-(-l)"' 



no-*;) 



x;(i-x ( o 



and x' is defined as: 



, J V .v, = 0 
■ • /- \l-.xj, Sj =\- 

As an example, consider the one dimensional case (i.e., TV = 1). For this case, the 
interpolation function is: 

yj [x] = u 0 [x] y[0] + 0 [xM0] + u } [*]j>[l] + vjx]a>;[l] 
and the interpolation basis functions are: 

u 0 [x] = (1 - 3x, 2 + 2x, 3 ) u } [ x] = (x, 2 ) (3 - 2x, ) 
\o M = 0 - *i ) 2 *i v 1?1 [x] = (-l)(x! ) 2 (1 - X, ) 

which gives: 

y, [x] = (1 - 3x, 2 + 2x, 3 ) j[0] + x, 2 (3 - 2x, ) j[l] 
+x ] (l-x 1 ) 2 a ) ;[0]-x l 2 (l-x 1 )ay[l] 



20 
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Similarly, for the two dimensional case (i.e., N = 2), the interpolation function is: 

J/[*l>*2] = 



v,,s 2 =0,l • 



or: 



«o,o[ x p x 2 1 y[Q, 0] + v 100 [x, ,x 2 ] d l y[0, 0] + v 2 0 0 [x, , x 2 ] 3 2 j;[0, 0] + 
"p,i t x i > x 2 1 ^[0, 1] + v 1>0t1 [x, , x 2 ] d, j[0, 1] + v 2 o , [x, , x 2 ] d 2 y[0, 1] + 

«1,0 [ X l » X 2 ] .PP. 0] + V l,I,o[ X 1 > X 2 3 d l yW\ + V 2,l,0 [*1 » X 2 ] 5 2 M 1 ' °] + 

[*i , x 2 ] Jil. 1] + v, , , [x, , x 2 ] a, y[\, 1] + v ? , , [x, , x 2 ] 5 2 y[\, 1] 
and the interpolation basis functions are: 

M)| .Jx„x 2 ]-(l- x[)(l x;)(l + x,'(l-2xr)x:(l 2x 2 )) 

and: 



which expand to: 



and: 



v, v , [x, ,x 2 ]= (-\f (i - *o (i - x; ) x; (i-x;> 

"o,o t x i > x i 3 = 0 ^ x i ) 0 - * 2 ) 0 + x i 0 - 2 x, ) x 2 (1 - 2 x 2 )) 

W 0.l[ X l' X 2] = ( 1 - X l) X 2( 1 + X l( 1 - 2x i)( 1 - X 2)( 2x 2- 1 )) 

M, 0 [x,,x 2 ] = x, (l-x 2 )(l + (l-x,)(2x l -l)x 2 (l-2x 2 )) 
[^i , ^ 2 ] = x, x 2 ( 1 H- (1 - x, ) (2 a:, -1)(1-x 2 )(2x 2 -1)) 

V l,0,0 [ X l > X 2 ] = 0 " X 1 ) 0 - X 2 ) X l 0 - X l ) 
V, p,l [ X ] » X 2 1 = 0 ~ X l ) X 2 X l 0 - X l ) 

, 10 [x p x 2 ]=-x 1 (l-x 2 )(l-x,)x l 



15 



^1,1,1 [^1 J "^2-1 ~ X2 (1 ) Xj 
V 2,0,0 t X P X 2 1 = 0 - *1 ) 0 ~ X 2 ) X 2 C 1 * X 2 ) 

V 2,l ,0 [*1 7 X 2 ] = X l 0 - *2 ) *2 (1 ~ *2 ) 
V 2,l,l [^1 ' *2 ] ~ ~ X ] X 2 0 — ^2 ) X 2 

As described previously, the reduced multicubic interpolation is typically used to 
interpolate a function FfX] where X (in the scatterometry case) is typically a vector or 
matrix that includes computed reflectance quantities such as complex reflectance 
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corresponding to multiple wavelengths, incidence angles, and polarization states. As 
shown in Figure 3, interpolation of X begins by locating an interpolation cell m in which 
X \s included (see step 302). The interpolation cell m is N-dimensional and includes X in 

each dimension (i.e., x] mj } < X f < X)'" j +1) ; m } = 0 . . . M } - 1; j = 1 . . . N ). 

5 In step 304, the interpolation argument X is mapped into a transformed argument 

jc. The transformation is equivalent to moving the interpolation cell m to a unit cell with 
its base vertex at the origin. The transformed X argument, denoted x = (x } ,x 2 ,...x N ) 9 is 
defined by: 

X = J - I • / =r 1 N 

X j " ^(jh/+1) ) ' 

1 0 Once X has been transformed, the problem of interpolating Y becomes a problem 

of interpolating a function y within the unit cell. In step 306, the values for y at each 
vertex of the unit cell are retrieved from the database. For the one dimensional case 
described above, this means that y[0] and j[l] are retrieved. For the two dimensional case, 
four y values would be retrieved: j[0,0], ^[0,1], >>[1,0] andj?[l,l]. 

15 In step 308, first order derivatives are obtained for y at each vertex. For the one 

dimensional case described above, this means obtaining dy[0] and dy[l] . Some 
implementations store the first derivatives as part of the interpolation database. In these 
cases, obtaining the first order derivatives is performed by database lookup. In other 
implementations, the first order derivatives are not stored and must be generated on-the- 

20 fly. Typically, this will involve a numerical method for estimating derivatives using the j> 
at each vertex of the unit cell as well as nearby y values (e.g., y[-l] and y[2] for the one 
dimensional case). 

In step 3 1 0, the interpolation function is evaluated using the y values and their 
associated first order derivatives. 

25 

Reduced-Cubic Interpolation with Finite-Difference Derivatives 

In practice, the derivative terms d i y[s] in Equations 16 and 17 may not be 

readily available and must be evaluated by finite differencing. In these cases, the present 
invention provides a finite difference approach to reduced-cubic interpolation. For this 
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approach, the gird-sampled function values y[s] are used for derivative estimation. 
Assuming that the sampling grid is uniform in each dimension (i.e., x ( p +l) - X { j* j) is the 
same for all m } = 0. . M } - 1 ), a simple centered finite differencing approximation can be 
used: 

5 5,y[-v]~^0'[.vtc',]-jlv-e ( l) (19) 

wherein e i is the unit vector: 

IX i = J 
[0, i * j 

If grid cell m is at the boundary of the data sampling range (i.e., m i - 0 or 
m i = M i - 1 for some / = 1 . . . Af ), then one of the terms y[s ± e f ] in Equation 1 9 will be 
10 outside the sampling range. In this case a one-sided derivative estimator may be used: 

5, y[s] =± (-±y[s ± 2e,\ + 2y\s ± *,] - f tfs]) (21) 

However, Equation 21 is less accurate than Equation 19, so the data sampling 
range should preferably be "padded" with data points external to the interpolation range 
to avoid the need for one-sided derivative estimators. (The padded sampling range should 
1 5 be defined such that for each grid point X (m) in the interpolation range, all points 
X {m±€i) for i = l...N are sampled.) 

A higher-order derivative estimator may be used to improve interpolation 
accuracy: 

diyW^iyis + el-yls-eV-^iyis + lel-yis-left (22) 

20 For grid cells at or near the sampling boundary, higher-order, asymmetric 

derivative estimators such as the following may be employed (but at the expense of 
compromised accuracy), 

{+ly[s±e l ]-^y[s±2e l ] + ^y[s±3e l ] / 
d, y[s] = ±{^±y[s - 1 y[s] + y[s ±e,) - ±y[s +2^]) (24) 



(23) 
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( 25 



d i y[s] = ± 



(25) 



-3y\s±2e i \+*y[s±3e l \-\y\s±4e i \) 
d l y[s] = ±(-%y[s] + 3y[s±e i ]-±y[ S ±2e l ] + ±y[s+3e,}) (26) 

The need for these asymmetric derivative estimators can be avoided by padding 
the sampling range such that for each grid point X (m) in the interpolation range, all 
points X (m± ° and A r(m±2e ' ) for i-l...N are sampled. 



10 



20 



Reduced-Cubic Interpolation with Point Matching 

With this method, the function matching conditions are similar to Equations 1 6 
except that finite differences are substituted for the derivative terms: 

2 O', l> + I - .v, I v - <vl ) = t 01 s +<?,\- -vl * - I); 

sr. =0,1; i,j = \...N 

The interpolation function is: 



i:.v y =0,10=1. f)\ 



M...N 



The interpolation basis functions, u s [x] and s [x] , are defined by: 

f V \ 



",[*] = 



Ho-*;) i + i 2>;(i-2*;.) 

..N 

x]{2-x]) 



s 



15 and x' is defined as: 



(27) 



(28) 



(29) 



•> \\-x p Sj=V 

If grid cell m is at the boundary of the data sampling range, Equation 28 is not 
applicable because it uses y data points outside the sampling range. In this case, the 
interpolation function may be defined by simply extrapolating from the nearest interior 
grid cell. However, if more than one of the X. coordinates is at the sampling limit the 

extrapolated function may not actually fit all the boundary points; and in any case, 
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accuracy will be compromised in the boundary cells. Thus the data sampling range 
should preferably be padded to avoid the need for extrapolation. 

Implementation Approach 

5 The primary difficulty encountered in implementing the above algorithms is in 

translating the multidimensional array notation into conventional programming language 
syntax. For example, the data samples F[X (m) ] could be represented in Fortran by 
multidimensional array elements Ysample (m(l) , m(2) , ... m(N) ), wherein m {j) = m j . (If 

Y is a vector or matrix, the Ysample elements could be pointers to Y data structures.) 
10 However, this approach is very constraining because the number of array dimensions (n ) 
must be a compile-time constant and cannot exceed seven (due to an inherent Fortran 
limitation). In C++, multidimensional arrays are represented by nested pointer 
indirection, e.g., Ysample [m[0]] [ m[i] ]... [m[N-i] •], wherein m(j) = mj x . Butthis 

approach is also impractical because the number of dimensions is predetermined by the 
1 5 Ysample data type and the multiple pointer indirections are computationally inefficient. 
A more practical approach is to represent all multidimensional arrays 
computationally as one-dimensional arrays, using a Fortran-style memory layout, and 
perform the array indexing operations manually. For example, the datum Y[X {m) ] would 
be represented computationally as 

20 Y[X (m) ] = Ysample [m__index] (30) 

wherein m i ndex is a scalar array index, which is a linear function of m ; 

m index = 

offset + stride [0]*m[0] + stride [ 1]* m[l] (31) 
+ ... stride [N-l]*m[N-l] 

mlj) = m J+] (32) 

The program logic, as well as the data structures, can be adapted for variable N. 
25 For example, following is a multidimensional array initialization routine using a 
conventional C++ array representation, 
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for (m[N-l] - extent [N-l] ; --m[N-l] >= 0; ) 

for (m[l] = extent [1] ; --m[l] >~ 0; ) . (33) 

for (m[0] 5= extent [ 0 ] ; - -m [ 0] >= 0 ; ) 
Ysample[m[0] ] [m[l| J...[m[N-l] ] 

(Following Fortran nomenclature, the " extent " array represents the dimension 
extents.) This code structure is not applicable to variable n , but with the Equation 30 
representation the multidimensional iteration can be implemented as follows for variable 

5 . N, 

bool done = false; 
int m__index = offset; 
for (int i - 0; i < N; ++i) { 
if (extent [i] == 0) { 

done ^ true; 

break; 

- } 

m[i] =? extent [i]-l ; 
in_index += stride [ i ]* m [ i ] ; 

} ^ . 
while ( (done ) { 

Ysample [m_index] = ...; 

done = true; 

for (int i - 0; i < N; ++i) { 
if (--m[i] >= 0) { 

ni_index - = stride ti];. 

done = false; 

break; 

} 

m[i] == extent [i] - 1; 
m_index + stride [i]* m[i] ; 
} " . .. 

1 (34) 



Comparative Computational Performance 

In each of the interpolation algorithms the interpolation function y,{x) is 
10 represented as a linear combination of the grid-sampled y values and (possibly) 

derivatives d i y , wherein the linear coefficients are the x -dependent interpolation basis 
functions. Typically, y is a large data structure and the computational overhead 
associated with calculating the basis functions (Equation 10, 14, 18, or 29) is dominated 
by the computation of the linear combination (Equation 9, 13, 17, or 28). In this case the 
1 5 computation runtime is approximately proportional to the number of summation terms in 



TWI-32310 



- 23- 



PATENT 



the linear combination. For multilinear interpolation (Equation 9), the number of 
summation terms is 2 N . For full multicubic interpolation (Equation 14) the number is 
A N , and for reduced-cubic interpolation (Equation 17) it is (1 + N)2 N . Thus, for large 
N the reduced-cubic interpolation has a very substantial efficiency advantage over full 
5 multicubic interpolation. 

The reduced-cubic interpolation algorithm with point matching (Equation 28) has 
(1 + 2N)2 N summation terms, as written. However, there are only (1 + N)2 N distinct y 
sample values appearing in Equation 28, so in practice this equation can be reformulated 
so that it also has (1 + N)2 N terms. This same condition also applies to Equation 17 with 
10 substitution of finite-difference derivatives, Equation 1 9. If the higher-order derivative 
estimator, Equation 22, is used in Equation 17, then the number of y samples per 

interpolation is (1 + 2N)2 N . 

Continuity and Smoothness 

15 Multilinear interpolation yields a global interpolation function Yj[X] that is 

continuous across grid cell boundaries, but which is not generally smooth (i.e., 
continuously differenti able). With full, multicubic interpolation, Y { [X] is globally 
continuous and smooth. The reduced-cubic interpolation function is globally continuous. 
It is not generally smooth across grid cell boundaries, although it is smooth at the grid 

20 points X (M) . The same condition holds when finite-difference derivatives are used, 

although in this case the derivatives d j Y l will not exactly match d. Y at the grid points. 

Reduced-cubic interpolation with point matching yields an interpolation function that is 
globally continuous but is not generally globally smooth (even at the grid points). 
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Interpolation Accuracy 

In characterizing interpolation accuracy, it is useful to consider a function y[x] 

represented by a Taylor series near x = 0 , 



( f \ 

(35) 

j=\...N 



yw= z n ji s / j m n 



5 wherein p = ( p } , p 2 , . . .p N ) and 0 = (0, 0, ... 0) . Each monomial term Y\j ( x jY J i n 
Equation 35 has an associated "monomial order" I p defined by . 

ZP=ZPj (36) 

7=1. ..N 

"Interpolation accuracy order" is defined as the largest integer n such that any 
monomial of order n or less has zero contribution to the interpolation error: 

10 yj[x]-y[x]. 

The accuracy order is related to the scaling relationship between the interpolation 
error and the grid sampling density. Suppose that the interpolation sampling interval 
sizes are all proportional to a common scaling factor s : 

Xp^-X^^eA^; m j =0... M } -1;- j = l...N (37). 
15 Then, based on Equations 7 and 37, each 3 y term in Equation 35 carries with it an 

implicit factor of e , and each monomial term (Xj) Pj has an implicit factor of s^ p 

multiplying it. Thus, the lowest power of e appearing in the interpolation error 
y } \x\^y\x\ (represented as a Taylor series in £-)willbe , wherein n is the 
interpolation accuracy order. 

20 Multilinear interpolation (Equation 9) has accuracy order one. Multicubic 

interpolation (Equation 13) and reduced-cubic interpolation (Equation 17) both have 
accuracy order three. If the finite-difference approximation of either Equation 19 or 21 is 
employed in Equation 17 the accuracy order is reduced to two; however any of the 
higher-order derivative estimators, Equations 22-26, preserves order-three accuracy. 

25 Reduced-cubic interpolation with point matching also has order-three accuracy. 
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Following are lowest^order approximations for the interpolation error, for each 
interpolation algorithm. 

Multilinear interpolation, Equations 9, 10: - 

;,M-jM = |&;#]^(i^) (37) 

Multicubic interpolation, Equations 13, 14: 

J>,M^*]£-i Sa, 4 jl« (x,) 2 (l-x,) 2 (38) 

i=\...N 

Reduced-cubic interpolation, Equations 17, 18: 

yi [x]-y[x] = ^ Zd.'yiO] (x,) 2 (l-x,) 2 



i=\...N 



i Z E^/jW *,o-*,)*,<j-*/> 

• , 7-2... ;V 



(39) 



Equations 17, 1 8 with the finite-difference derivative, Equation 19, substituted in 
10 Equation 17: 

^M-^[xj = l Z d ' 3 yW x ( .(l-x ( )(l-2x ( .) (40) 



i=\'.:.N 



Equations 17, 1 8 with a higher-order finite-difference derivative, Equation 22: 
same as Equation 40. 

Reduced-cubic interpolation with point matching, Equations 28, 29: 

y,[x] -y[x] = Z d ? 0 + 0 - x,)(2 - x/) 
15 i=l - JV (42) 

-i S E^a/^] x,(i-x,)x,o-x,) 

j=2...N /=t...y-l 

[Note: Equations 40 and 42 have a similar form; but the factor (x,) 2 (1 - x,.) 2 

multiplying d, 4 y[0] in Equation 40 takes on a maximum value of ^~ over the unit 

interval 0 < x, < 1 , whereas the corresponding factor (1 + x,)x, (1 - x,)(2 - x ; ) in 

Equation 42 takes on a maximum value of . Thus, the accuracy of the point-matching 

20 method is an order of magnitude worse than the method employing accurate derivatives.] 
The above error estimates do not consider the case where asymmetric derivative 
estimators (Equations 21, 23-26) are used. However, it is evident from the following 
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lowest-order estimates of the derivative estimators that a significant accuracy penalty is 
incurred by using asymmetric estimators: 

Equation (19): 

5 13 (43) 

= id, 3 y[s] 



Equation (21): 



±(- \ y[s ± 2e, ] + 2y[s ± <?,]-§ y[s]) - d, y[s] 



(44) 



Equation (22): 



(46) 



^(y[s + e l ]-y[ S -e i ])-^(y[ S + 2e l ]-y[s-2e i ])-d l y[s] 

(45) 

10 Equation (23): 

±(-ly[s + e l ]-5y[ S ]+ly[s±e i ]-\y[s±2e l ] + J 2 y[ S ±3e i ]yd,y[s] 
~^y[s\ 

Equation (24): 

±(-\y[ S + e i yy[s] + y[ S ±e I ]-^y[ S ±2e l ])-d i y[s} 

Equation (25): 



=*M 4 y[s] 



(47) 



i5 ±(-^y[s] + 4y[ S ±e l ]-3y[s±2e l } + ^y[ S ±3e l ]~ly[s±4e l ])~d l y[s] 

- - I <V M 



Equation (26): 

±(- l iy[s] + 3y[s±e l )-ly[ S ±2e i ]+ly[ S +3e l ])-d i y[ S ] 



(49) 
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Discussion 

The multilinear and multicubic interpolation functions (Equations 9 and 13) are 
obtained by applying one-dimensional interpolation in each dimension separately. For 
example, applying linear interpolation successively to x, , x 2 , . . . , yields a multilinear 
5 interpolation function of the form 

p:pj*QMj=?l-N) j=\...N ' 

wherein the 2 N coefficients a p are determined from the corresponding 2 N point- 
matching conditions, Equation 8. Substituting yj[x] from Equation 50 into Equation 8, 
yields a set of 2 N equations that can be solved for the a p coefficients. a p will evidently 
10 be a linear combination of the constraint values y[s] appearing in Equation 8; hence 
substitution of the solution in Equation 50 yields a expression for y } [x] as a linear 
combination of the y[s] constraints, as expressed in Equation 9; wherein the u s [x] terms 

are multilinear functions of x comprising the same monomial terms appearing in 
Equation 50. 

15 Rather than solving for the a coefficients explicitly, the basis functions u s [x] 

are obtained by direct substitution of Equation 9 into Equation 8 (with jc = s'), 

Z " s [s'M^y[s'l> ^=0,1; j = l'...N (51) 

j;A V =0,iO=l...N) 

This applies to any function y ; hence the basis functions satisfy the condition 

* = V ^=<U 7 = 1.. .N (52) 
[0 s*s 7 

20 W 5 IXI is a 2 N -term multilinear function (similar to Equation 50) whose 

coefficients are determined by the 2 N conditions defined by Equation 52. (There are 2 N 
distinct s' values in Equation 52.) Based on the symmetry of Equation 52, the following 
condition holds for any j = 1 . . . N : . 

»-:*!, "H,, , .v < 53) 
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Thus, it suffices to consider the case s = 0 . Given the basis function u 0 [x] , Equation 53 
can be applied to obtain u s [x] for any s ( ^ f = 0, 1; j = 1 . . . N ). Assuming multiplicative 
separability of in x, , x 2 , ... x N , the following result is readily obtained from 

Equation 52: 

5 u 0 [x]= llC -*,) (54) 

j??\...N 

and application of Equation 53 yields the result, Equation 10. 

The methodology outlined above can also be applied to multicubic interpolation. 
In this case, the interpolation function has a form similar to Equation 50, except that p } 

ranges from 0 to 3: 

10 y[x] = y,[x]= £ "/I! (*/)"' < 55 > 

p:p y =0...3O=l : ..AT) .j=]...N . 

This expression, together with the constraints defined by Equation 12, implies an 
interpolation function having the form of Equation 13 wherein the basis functions u h s [x] 

comprise the same monomials appearing in Equation 54. The defining conditions for 

"*,,M are 

*} r '1 - J 1 ' k ' = kands ' = s 
J J UkASi '{0, AVAor*'*.v ' (56) 

k' J ,k J ,s' j ,s J = 0,\; j = \...N 
The basis functions satisfy the symmetry conditions, 



J 



Hence, given the basis functions u k 0 [x] , all the other basis functions can be determined 
by application of Equation 57. Assuming multiplicative separability of u k 0 [x] in x, , x 2 , 
20 ... x N , yields: 

wherein / is defined by Equation 15; and application of Equation 57 then yields 
Equation 14. 
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10 



15 



The starting point for reduced-cubic interpolation is the constraint conditions, 
Equations 16, which involve first derivative matching but do not involve the mixed first 
derivatives appearing in Equation 12. Since equation 1 6 defines (1 + ^)2* constraints, 
an interpolation function is needed that has this same number of free parameters. The 
following functional form meets this criterion: 



j=\...N )j=}...N 



(59) 



Equations 16 and 59 imply that y } [x] can be represented as a linear combination of the 
constraint values y[s] and d i y[s] , as represented in Equation 17; wherein the basis 
functions u s [x] and v, s [x] comprise the same monomials that appear in Equation 59. 

After back-substituting Equation 17 in Equation 16 and matching corresponding 
constraint coefficients, the following defining conditions for u s [x] and v, s [x] are 

obtained, 

, fl, s' = S 

[0, s *s 

d r u s [s'] = 0 (60) 
[l, /' = i and s' = s 
[0, i *iors'*s 
s'j, Sj =0,1; i\i,j = \...N 



The basis functions satisfy the following symmetry relationships, 



v t si x i 1 



l = J 

The following functions satisfy Equation 60 with s = 0 , 



(61) 
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u 0 [x] = 



( Y 

( \ (62) 



and applying Equations 61, 18 are obtained. 

Global continuity of Y f [X) can be demonstrated for multilinear interpolation as 
follows: If Equation 9 is evaluated on one of the interpolation cell's hyperplane 
5 boundaries, e.g. at x, = 0 , then the computed value y f [x] will depend only on y[s] 
values sampled in the hyperplane (i.e., with s x - 0). Hence the untransformed 

interpolation function T/ w) [ X] , evaluated on the hyperplane X x == X^" h) (on the 
boundary of cell m ), will depend only on Y[X] values sampled in the hyperplane. 
Another interpolation function Y^' l) [X] defined on adjacent cell m -e } , will similarly 

10 have a functional form in the hyperplane X } = X^ ih) that is defined entirely by the same 
Y[X] values; hence the two functions Yj {m) [X] and F / (m " e,) [X] will be matched at the 
common boundary between cell m and cell m-e } . The same condition applies to all 
bounding hyperplanes between interpolation cells; hence the global interpolation function 
Yj[X] is everywhere continuous. 

15 The above line of reasoning can be extended to show that all of the interpolation 

algorithms outlined in Section 3 yield globally continuous interpolation functions, and 
that multicubic interpolation yields an interpolation function that is also globally smooth. 
The reduced-cubic interpolation function, Equation 1 7, is not globally smooth; but it is 
smooth at the grid points because the function derivatives are constrained by Equation 1 6 

20 to be matched between adjacent grid cells. 

Assuming that no approximations are made in the interpolation constraints, the 
interpolation error >>/[*]- j[jc] will be identically zero (by construction) when y[x] is 

any polynomial comprising the same monomial terms that are in y { [x] . Thus, to 
determine an error estimate, it is appropriate to limit consideration to the lowest-order 
25 monomial terms in y[x] that are not in yj[x] . For example, the multilinear function 

yj[x] defined by Equation 50 includes all monomials up to order 1; thus the multilinear 
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interpolation error is dominated by monomials of order 2 that are not in the multilinear 
expansion. There are N such monomials: (x y ) 2 , j = 1. . . N . Considering, for example, 
the monomial function 

y[x] = ( X] ) 2 (63a) 

5 the terms y[s] in Equation 9 are 

fO, .v, = 0 

y[s] - j , . (63b) 

With this substitution, Equation 9 reduces to y, [x] = x, , 

y,[*]*= Z " S M = E ■ ■ S "(i,.v J ,...v ;V )W 

i:i,=l;.t y «0,l(7=2...yV) i 2 =( >>l Sjv-0,1 

= i ••• i x > n (i-^)=^ n (630 

" s 2 =0,1. .v, v =0,l J-2...N . ' j-2..,N .v 7 -0,l 

7 = 2. ..JV 

The interpolation error, for defined by Equation 63 a, is 

10 J/M-jM-x.O-^) (63d) 

Generalizing this to an arbitrary function y[x] defined by a Taylor series (Equation 35), 
the lowest-order interpolation error for multilinear interpolation is given by Equation 38. 

A similar procedure yields the interpolation error estimates for multicubic 
interpolation (Equation 39), reduced-cubic interpolation with exact derivatives (Equation 

1 5 40), and reduced-cubic interpolation with point matching (Equation 42). Each of these 
methods has the property that the interpolation error is zero for any monomial included in 
Vi [ x ] • F° r example, in the case of multicubic interpolation, the interpolation function 
yj[x] defined by Equation 55 includes all monomials of order 3 or less, and the only 
order-4 monomials not included in the multicubic expansion are (x y ) 4 , j = 1. . .N . Thus 

20 the multicubic interpolation error, Equation 39, can be estimated by considering the 
special case y[x] = (x^ 4 (analogous to Equation 63a). The reduced-cubic expansion, 
Equation 59, also includes all monomials of order 3 or less, but the order-4 monomials 
not included in Equation 59 include the functions (x y ) 4 and also (x,) 2 (x,) 2 ; i * j . 
Functions y[x] of both forms must be considered in deriving Equations 40 and 42. 
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Equation 41 is a special case because the interpolation algorithm relies on an . 
additional approximation, the derivative estimate (Equation 19), which is only accurate to 
order 2. (In this case, the y l derivatives are only constrained to match finite-difference y 
derivatives at the grid points.) The lowest-order monomials for which Equation 1 9 is in 
5 error are: (Xj)\ J = 1 . . . N ; thus the interpolation error estimate in Equation 4 1 is 

determined by considering functions y[x] of this form. However, if the higher-order 
derivative estimator, Equation 22, is used, the accuracy estimate in Equation 40 is 
unchanged because Equation 22 has accuracy order 4 (see Equation 45). 

10 Comparison 

The primary advantages and disadvantages of the interpolation algorithms 
considered above are summarized as follows: 

The prior art multilinear interpolation (Equation 9) is the most computationally 
efficient algorithm, requiring 2 N summation terms per interpolation. The interpolation 
15 function is globally continuous, but is generally not smooth at the grid cell boundaries. 
Multilinear interpolation only has accuracy order 1 (Equation 38). 

Full, multicubic interpolation (Equation 13) (also prior art) has accuracy order 3 
(Equation 39). (This assumes accurate derivatives are used.) The multicubic interpolation 
function is globally continuous and smooth. However, the algorithm is very inefficient, 
20 requiring 4^ summation terms. 

Reduced-cubic interpolation with exact derivatives (Equation 17) also has 
accuracy order 3 (Equation 40), but requires only (1 + N)2 N summation terms. The 
interpolation function is globally continuous. It is not generally smooth, except at the 
grid points. A disadvantage of this method (as well as the multicubic interpolation 
25 method considered here) is that it requires derivative information. 

The need for exact derivatives can be eliminated by using finite-difference 
derivatives, but there are tradeoffs. If a simple, centered-difference derivative estimator 
(Equation 19) is used, the interpolation accuracy is reduced to order 2 (Equation 41). A 
higher-order derivative estimator (Equation 22) may be used to restore order-3 accuracy, 
30 but then the number of summation terms is increased to: (1 + 2N)2 N . 
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The above tradeoff can be circumvented by using the reduced-cubic interpolation 
method with point matching (Equation 28). This method does not require derivatives; it 
has accuracy order 3 (Equation 42); and the required number of summation terms is 
(1 + N)2 N . However, the accuracy error is typically much larger (e.g., by an order of 
magnitude) than that of the derivative-matching methods (assuming that accurate 
derivatives are used). Another possible disadvantage is that the interpolation function 
may have poorer smoothness characteristics compared to the other reduced-cubic 
methods. (The point-matching interpolation function is generally not smooth between 
grid cells, even at the grid points.) 
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